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Abstract 

The unsteady structure of a supersonic jet is highly 
three dimensional, though the mean flow is axisymmet- 
ric. In simulating a circular jet, the centerline represents 
a computational boundary. As such, spurious modes can 
be generated near centerline, unless special attention is 
given to the behavior of the 3D structure near the center- 
line. Improper treatment of the dependent variables near 
the centerline results in the solution diverging or being 
suitable only for small amplitude excitation. With a care- 
ful treatment of the centerline formulation, no spurious 
mode is generated. The results show that a near linear 
disturbance growth is obtained, as die linear stability the- 
ory indicates. At high levels of excitation, nonlinear de- 
velopment of disturbances is evident and saturation is 
reached downstream. 

1. Introduction 

Jet noise suppression has become a critical issue for 
the development of high speed civil transport plane. The 
jet noise is generated by the time dependent flow fluctu- 
ations in the near field which are associated with pres- 
sure fluctuations that propagate to the far field producing 
the radiated sound. Experiments have shown that the 
measured sound fields appear to emanate from a region 
about 10 diameters downstream of the nozzle exit [1]. 
This noise-producing initial region of the jet is character- 
ized by a large scale vortical structure and can be viewed 
as having a wavelike nature. It is believed that the large 
scale structure is more efficient than the small scale 
structure in radiating sound [1-5]. This indicates that die 
initial development of the jet should be clearly resolved 


* Senior research associate. Member AIAA 

** Senior scientist and technical leader, CAA, Associate 

fellow AIAA 


so that an accurate noise prediction can be made. The 
near flow field is described by the unsteady Navier- 
Stokes equations. However, direct numerical simulation 
can not resolve all scales of motion for high Reynolds 
number flows. It is a pp r opri ate to perform large-eddy 
simulations to accurately capture the large scales of mo- 
tion while modelling the sub-grid scale turbulence. 

The use of large-eddy simulations (LES) as a tool 
for prediction of the jet noise source has been proposed 
by Mankbadi et al. [6], Not only the mean flow must be 
calculated accurately, but also the physical flow fluctua- 
tions must be accurately predicted since the sound source 
is given in terms of the flow fluctuations. Since the com- 
putational domain is usually finite, the numerical bound- 
ary conditions can generate spurious modes that render 
the computed flow fluctuations totally unacceptable. In 
simulating a circular jet, the centerline (r=0) represents a 
computational boundary. Boundary condition for axi- 
symmetric mean flow is obvious. However, the unsteady 
structure of a supersonic jet is highly three dimensional, 
even though the mean flow is axisymmetric (see, for in- 
stance, Mankbadi [5], Michalke [7]). As such, spurious 
modes can be generated near the centerline, unless spe- 
cial attention is given to the behavior of the three dimen- 
sional structure near the centerline, which is the subject 
of the present wok. 

The numerical solution of the Navier-Stokes equa- 
tions in cylindrical coordinates requires the proper treat- 
ment of the discretized equations at the centerline. In 
axisymmetric jet flow, Mankbadi et al [6] derived a new 
set of equations at the centerline from the original equa- 
tions using L’Hospitals rule to circumvent numerical 
problems associated with the geometric singularity in the 
formulation. In the present study, three approaches to the 
centerline treatment are considered, namely, asymptotic, 
averaging and interior points approaches. The effect of 
sub-grid scale turbulence stresses is not taken into ac- 
count in this study to avoid any uncertainty from the tur- 
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bulence models. A supersonic jet with Mach number = 

1_5, Reynolds number = 1.27 x 10^ based on the nozzle 
exit diameter and jet centerline parameters is considered. 
The outer stream is 0.25 of the jet exit velocity, and the 
jet temperature ratio is 0.5. Time -harmonic disturbances 
are imposed at the inflow boundary of the jet, and the 
subsequent development of disturbances are examined. 

2. Governing Equations 


The flow field of a supersonic jet is governed by the 
compressible Navier-Stokes equations, which can be 
written in cylindrical coordinates as 
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Here Q is the unknown vector, F, G, and H are the fluxes 
in the x, r, and 4 directions, respectively; S is the source 
term that arises in cylindrical polar coordinates; and k is 
thermal conductivity. The total enthalpy is I, the total en- 
ergy is E, and a^ are the viscous stresses. This system of 
equations is coupled with the equation of state for a per- 
fect gas. 


3. Numerical Scheme 


where 

Q = [p, pu, pv, pw, pE] T 
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The importance of the dispersion and dissipation of 

(2) a given scheme, which used in connection with compu- 
tational aeroacoustics, were highlighted by Hardin [8]. 
Both effects are crucial in computational aeroacoustics, 
and can render the computed unsteady part of the solu- 
tion completely unacceptable. As such, high -order accu- 
rate schemes are required for problems in computational 
aeroacoustics. 

(3) 

A fourth -order accurate in space, second-order accu- 
rate in time scheme is used, which is an extension of the 
MacCormack scheme by Gottlieb and Turkel [9]. Mank- 
badi et aL [6] used this scheme to study the structure of 
axisymmetric supersonic jet flow and its radiated sound. 
Ragab and Sheen [10], and Farouk, Oran and Kailasan- 
ath [11] have also successfully applied this scheme for 
the study of nonlinear instability problems in plane shear 
layers. Sankar, Reddy and Hariharan [12] performed a 
comparative study of various numerical schemes for 

(4) 

aeroacoustics applications, and found that this scheme 
offers high spatial accuracy. In this scheme, the operator 
is split into three one-dimensional operators and applied 
in a symmetric way to avoid biasing of the solution: 

QD + 2 = L 2x L 2r L 2$ L l$ L lr L lx QD ^ 

where L represents the one-dimensional operator. Each 
operator consists of a predictor and a corrector steps, and 

( 5 ) each step uses one-sided differencing: 
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and likewise for the radial and azimuthal directions. The 
scheme becomes fourth-order accurate in the spatial de- 
rivatives when alternated with symmetrical variants. Let 
L] be the one-dimensional operator with forward differ- 
ence in the predictor and backward difference in the ax- 
rector, then Lj will be the one-dimensional operator with 
backward difference in the predictor and forward differ- 
ence in the corrector. 


4. Boundary Conditions 


The scheme uses one-sided differences for the flux- 
es. A cubic extrapolation is used to obtain the fluxes at 
two ghost points outside the computational domain in or- 
der to update the boundary points, i.e. 

F m-l= 4F m- 6F n,-l + 4F m-2- F m-3 < 10 > 

F m + 2 = 4F !D + l -6F m + 4F m-l -F m-2 (11) 

The physical outflow boundary conditions for the com- 
putation are derived using linearized characteristics 
[13,14] to permit the unsteady flow properties to pass 
without producing non-physical reflections. A summary 
of the outflow boundary conditions is given in Mankbadi 
etal. [6]. 


The inflow was taken to be the mean flow with a hy- 
perbolic-tangent profile plus time-harmonic disturbanc- 
es. The initial mean axial velocity profile is given by 
Michalke and Hermann [IS] as 



( 12 ) 


where uj is the exit velocity of the jet, Uo denotes the ve- 
locity of coflow and 0 is the momentum boundary layer 
thickness of the jet shear layer. The momentum thickness 
6 and radius r are normalized by the nozzle exit radius R. 
The corresponding temperature is specified by the Croc- 
co’s relation, which can be expressed as 


where M is the Mach number. The time-harmonic distur- 
bances consist of a single, symmetrical pair of the first 
helical mode, i.e. 

[u* v’ w’ p' p] = «t»(r)e l(aX_<0t) cos(^) +CQ14) 

where a is the wavenumber and CC is the complex con- 
jugate. The angular frequency co is chosen such that the 
corresponding Strouhal number Sfc=f(2R)/uj=0.125. The 
eigenfunctions <l>(r) are obtained from the solution of lin- 
ear stability equations for velocities, density, and pres- 
sure respectively. 

^(0 = [u(r) v (r) w(r) p(r) p(r)] ( 15 ) 

We now focus our attention on the flow behavior near 
the centerline, as r approaches zero. 

S. Centerline Formulations 

In implementing the fourth -order scheme near the 
centerline, two ghost points will be needed for updating. 
We found that using two ghost points near the centerline 
generates nonphysical oscillations and causes the code to 
blow up. Instead, we used the second-order MacCor- 
mack scheme for updating this point Therefore, only one 
ghost point is needed. Three approaches to centerline 
treatment are considered and discussed below. 

S.l Asymptotic Behavio r of Navier-Stokes Equations 
Near the Centerline 

Starting from the Navier-Stokes equations, and re- 
quiring that the flow variables should be finite for r=0, 
we find that the azimuthal velocity w must satisfy the dif- 
ferential equation 

2 

*^ + w= 0 (16) 

a * 2 

The solution for w near the centerline can thus be written 
as 

w = Acos$ + Bsin$ (17) 

We chose to study a top-bottom symmetrical jet, Le. the 
three-dimensional helical modes come in pairs as the ex- 
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perimental investigation of Cohen and Wygnanski [16], 
and the theoretical investigation of Mankbadi [5,16] 
have shown. In this case, the constant A in the above 
equation is set to zero. The formulations for u, v, p and p 
can be derived in the same way. It follows from the Navi- 
er- Stokes equations that the flow behavior at r=0 can be 
written as 

wW * 

*«> - 

^ < is) 

£*♦> - *<♦.)=?„ 

These equations are valid at r=r 0 where r 0 is a small 
quantity and is taken to be 0.01 in the present study. Thus 
the above equations describe the solution at 6 on the cir- 
cle r=r 0 in terms of the solution at 4> 0 . In implementing 
this condition, the solution at is obtained by taking the 

average of the values at 0° and 180°. 

S2 Averaging 

Another alternative that can be used is the simple averag- 
ing. The idea is that the flow variable at centerline r=0 is 
single-valued and independent of the azimuthal direc- 
tion. In implementing this approach, the flow variables u, 
v, w.p, and p at r=0 are taken to be the azimuthal average 
of all points at r=Ar. 

S3 Centerline as an interior point 

The centerline boundary condition is in fact an arti- 
ficial one that arises from the fact that polar coordinates 
are used with r starting from zero. Imagine performing 
full jet simulation with r extending from — to °° , thus 
point “a” closest to the centerline (figure 1) can be treat- 
ed as an interior point. If a second-order scheme is used, 
then the flux at point "b” is needed for updating point 
“a”. Since we are solving half a jet, point “b” is a ghost 
point. But because of symmetry, the flux at point **b" can 
be related to that of the interior point “c” as follows 

M. = M (19) 

b c 

N. = -N (20) 

b c 


where M represents the mass, axial or radial momentum, 
or energy flux, while N is the momentum flux in the az- 
imuthal direction. 

6. Results and Pise usrion 

A supersonic jet with Mach number 13, based on 
the jet exit parameters, is considered in this study. The 
velocity ratio UoAij is 0.25, and the temperature T</Tj is 

0.5. The Reynolds number based on the jet diameter and 

jet exit centerline velocity is 1.27X10 6 . The parameter, 
momentum thickness 6, in the initial axial velocity pro- 
file is 0.125. The inflow is taken to be the mean plus the 
time harmonic disturbance. 

Q = Q + e • Q’ (21) 

where Q and Q’ were given in the previous section, and 
e is the input excitation level. The computational domain 
(as shown in figure 2) extends 5 radii in the radial direc- 
tion, and 50 radii in the axial direction. 6 ranges from 0 
to 180 degrees in the azimuthal direction, since we con- 
sider a top-bottom symmetrical jet, and only half of the 
circular jet is solved. The computational grid consists of 
300x100x13 mesh points in the axial, radial and azi- 
muthal directions respectively. The mesh is uniform in 
the axial and azimuthal directions, and stretched in the 
radial direction with concentration of the grid points 
around r=R. The computations started with a time step 
corresponding to a CFL number of 0.3, and later on, after 
the initial transient purged out of the computational do- 
main, a finer time step was used to collect the computed 
data for use in Fourier transform. 

Figure 3 shows a snap shot of the predicted density 
contours using the asymptotic formulation, and a clean 
behavior near the centerline is evident. The excitation 

level £ is lxl0‘ 5 in this case. For the sake of comparison, 
figure 4 shows a snap shot of the density with no careful 
treatment of the centerline. These results were obtained 
using a second-order MacCormack scheme for the points 
at r=r 0 +Ar, and second-order extrapolation for the points 
at t=To. As indicated in figure 4, nonphysical values were 
generated near the centerline. Figure 5 shows the spectra 
at r=r 0 , x=40R. This figure shows that the amplitude 
peaks at the forcing frequency St=0.125 and its harmon- 
ic, and no spurious modes were generated when the as- 
ymptotic form ulatioo is used. Without careful treatment 
of the centerline, spurious modes were produced that 
rendered the solution totally unacceptable. 
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Figure 6 shows the growth of the input axial velocity 
disturbance at r/R=l for several excitation levels using 


the asymptotic formulation. The time span in which the 
data was collected is 16 characteristic times, where the 
characteristic time is defined as the ratio of nozzle exit 
radius to the jet centerline velocity. This time span corre- 
sponds to the period of the input disturbance when the 
Strouhal number is 0.125. Fourier transforms are per- 
formed for the collected time dependent data to obtain 
the amplitude of the disturbances at each location. As 
one can see in figure 6, the disturbance grows linearly, as 
linear stability theory indicates. At high levels of excita- 
tion, nonlinear effects occur and saturation is reached 
downstream. 

Next, we show the results obtained using simple av- 
eraging at die jet centerline. A snap shot of the predicted 
density contours is shown in figure 7. Oean behavior 
near die centerline is obtained. Figure 8 presents the 
growth of the input axial velocity disturbances at r/R=l 
for various excitation levels. An initial linear develop- 
ment of disturbances can be clearly seen and nonlinear 
effects take place when moving further downstream. The 
saturation of disturbance growth is also reached down- 
stream. 

Figure 9 shows a snap shot of the predicted density 
contours for e=lxl0' 5 using the interior point approach. 
No nonphysical value is generated near the centerline ex- 
cept the region near the outflow boundary. Figure 10 pre- 
sents the growth of the input axial velocity disturbances 
at r/R=l for various excitation levels. The same phenom- 
enon as those of the previous two approaches occurred, 

i.e. an initial linear growth of disturbance followed by its 
nonlinear development The saturation of disturbance 
growth is also reached downstream. Figure 1 1 shows die 
comparison of the growth of axial velocity disturbances 
at r/R= 1 using the three different approaches. The results 
obtained using the interior point approach have a slightly 
higher growth rate than the results obtained using the as- 
ymptotic and averaging approaches. 

The three-dimensional structure of the flow field is 
shown in figures 12-14. Figure 12 shows the snap shot of 
the computed velocity vectors in r-$ plane at x=283, 
333 and 38.3. As demonstrated in this figure, the cross 
flow velocity increases and changes its direction when 
the flow goes downstream; a clear evidence of the devel- 
opment of stream wise vorticities. Figure 13 shows the 
kinetic energy and vorticity contours at $=0 and 180 de- 
grees. It is seen that the symmetry of the flow no longer 
exists downstream. Figure 14 shows the iso-surfaces for 
the kinetic energy and vorticity magnitude at time t=100. 
The helical nature of die structure and the roll-up of vor- 
tices are evident 
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Direct simulations of a supersonic round jet flow 
field were presented with emphasis on the numerical 
treatment of the centerline (r=0). Three approaches, as- 
ymptotic. averaging and interior point were considered 
in the present study. Similar results were produced by all 
three approaches; a clean behavior of the flow near cen- 
terline, linear growth followed by nonlinear develop- 
ment of the disturbances, and the helical natural of the 
flow structure. Currently, attempts are under way to 
patch the computed nonlinear flow field to a linearized 
Euler solution to obtain the far-field sound. 
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Fig. 1 Schematic of meshes near centerline, r*4> plane. 


Nonreflecting boundary condition 



Fig. 2 Computational Domain 
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(a) (b) 

Fig. 3 Snap shot of the predicted density contours at t=100 using asymptotic formulation, (a) r-<}> 
plane, x=40 (b) x-r plane, 4>=0 deg. 



Fig.4 Snap shot of the predicted density contours at t=100 using extrapolation, (a) r-<}> plane, x=40 
(b) x-r plane, <J>=0 deg. 




Fig. 5 Comparison of axial velocity spectra at Fig. 6 The growth of axial velocity disturbance 
the centerline, x=40, <f>=0 deg. at r=l using asymptotic formulation. 
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Fig. 7 Snap shot of the predicted density contours 
at t=100 using averaging approach, x-r plane, <j>=0 
deg. 
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Fig. 9 Snap shot of the predicted density contours 
at t=100 using interior point approach, x-r plane, 
<j>=0 deg. 



Fig. 8 The growth of axial velocity distur- 
bance at r=l using averaging approach. 
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Fig. 10 The growth of axial velocity distur- 
bance at r=l using interior point approach. 



Fig. 1 1 Comparison of the growth of axial 
velocity disturbance at r=l. 
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(c) 


Fig. 12 Snap shot of the computed velocity vectors at t=100 in r-<}> plane, (a) x=28.3 (b) x=33.3 (c) 
x=38.3. 
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Fig. 13 Snap shot of the predicted kinetic energy and vorticity contours at 0=0 and 0=180 deg., 
t=100, (a) kinetic energy, (b) vorticity magnitude. 



Fig. 14 Snap shot of the predicted kinetic energy and vorticity isosurfaces at t=100, (a) kinetic 
energy=0.5, (b) vorticity magnitude=0.8. 
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